(A=Table[a[i,j], {i,2}, {j,2}]) // MatrixForm

(B=Table[b[i,j], {i,2}, {j,2}]) // MatrixForm

ArrayFlatten[Outer[Times, A,B]] // MatrixForm

Unprotect[NonCommutativeMultiply];
A_ ** B_ := ArrayFlatten[Outer[Times, A,B]];
A**B // MatrixForm

MatrixForm /@ ({Id, s1, s2, s3} = {(({{1, 0}, {0, 1}})), (1/2)(({{0, 1}, {1, 0}})), (1/2)(({{0, -I}, {I, 0}})), (1/2)(({{1, 0}, {0, -1}}))})

b[x_,y_] := x.y-y.x;
MatrixForm /@ {b[s1,s2], b[s2,s3], b[s3,s1]}

(H=4((1/4)Id**Id-s1**s1-s2**s2-s3**s3)) // MatrixForm

MatrixForm /@ {b[H, Id**s1+s1**Id], b[H, Id**s2+s2**Id], b[H, Id**s3+s3**Id]}

(H2=2H) // MatrixForm

Eigensystem[H2]

(H3=H**Id+Id**H+4((1/4)Id**Id**Id-s1**Id**s1-s2**Id**s2-s3**Id**s3)) // MatrixForm

Eigensystem[H3]

(H4=H**Id**Id+Id**H**Id+Id**Id**H+4((1/4)Id**Id**Id**Id-s1**Id**Id**s1-s2**Id**Id**s2-s3**Id**Id**s3)) // MatrixForm

MatrixForm /@ Eigensystem[H4]

Time this took: about 30 minutes.